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Scientific  Progress  and  Accomplishments 

contract /grant  number:  W911NF-09-1-0111 
ARO  proposal  number:  56031-EG 
December  31,  2012 
Richard  A.  Regueiro 
Associate  Professor 

Department  of  Civil,  Environmental,  and  Architectural  Engineering 
University  of  Colorado,  Boulder 

The  outline  of  this  section  of  the  final  report  is  as  follows: 

1.  overview  of  micromorphic  finite  element  (FE)  results  presented  in  various  published 
journal  articles  [Isbuga  and  Regueiro,  2011,  Regueiro  and  Isbuga,  2011,  Isbuga  and 
Regueiro,  2012]  and  PhD  thesis  [Isbuga,  2012]  funded  by  this  research  grant; 

2.  overview  of  discrete  element  and  micropolar  continuum  finite  element  coupling  in  one 
dimension,  and  numerical  results  [Regueiro  and  Yan,  2013]; 

3.  formulation  of  overlap  coupling  for  grain-FE  scale  and  micromorphic  continuum  FE 
scale  (two  methods); 

4.  overview  of  adaptivity  approach; 

5.  summary  of  future  research 

1  Overview  of  micromorphic  finite  element  (FE)  re¬ 
sults 

This  section  summarizes  briefly  the  micromorphic  finite  element  (FE)  results  presented  in 
various  published  journal  articles  [Isbuga  and  Regueiro,  2011,  Regueiro  and  Isbuga,  2011, 
Isbuga  and  Regueiro,  2012]  and  PhD  thesis  [Isbuga,  2012]  funded  by  this  research  grant. 
More  details  can  be  found  in  the  cited  articles  and  thesis,  provided  as  attachments. 

The  papers  Isbuga  and  Regueiro  [2011],  Regueiro  and  Isbuga  [2011]  focussed  on  the  finite 
strain,  materially  linear  isotropic  clastic,  three-dimensional  finite  element  implementation 
and  numerical  results  using  Tahoe.  We  showed  that  the  implementation  is  convergent  with 
respect  to  spatial  discretization  (i.e.,  mesh  size),  even  in  the  presence  of  a  boundary  layer. 
Upon  increasing  elastic  length  scale,  the  boundary  layer  was  found  to  become  less  sharp. 
For  pressure-sensitive  Drucker-Prager  elasto-plasticity  [Isbuga  and  Regueiro,  2012,  Isbuga, 
2012],  it  was  found  that  choice  of  micromorphic  elastic  parameters  in  conjunction  with  plastic 
strength  parameters  will  influence  the  stress  paths  and,  in  turn,  yielding.  The  boundary 
conditions  on  the  micro-displacement  tensor  3>/f  had  a  significant  influence  on  the  results.  A 
plastic  length  scale  term  has  not  yet  been  implemented,  and  it  is  anticipated  that  this  will 
influence  the  numerical  results  where  there  is  a  gradient  of  an  internal  state  variable  (ISV) 
in  a  region  of  large  deformation. 
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2  Overview  of  discrete  element  and  micropolar  contin¬ 
uum  finite  element  coupling  in  one  dimension,  and 
numerical  results 

The  ultimate  problem  of  interest  is  shown  in  Fig.l.  The  discrete  element  (DE)  and/or  finite 
element  (FE)  representation  of  the  particulate  micro-structure  is  intentionally  not  shown 
in  order  not  to  clutter  the  drawing  of  the  micro- structure.  The  grains  (binder  matrix  not 
shown)  of  the  micro-structure  are  ‘meshed’  using  DEs  and/or  FEs  with  cohesive  surface  el¬ 
ements  (CSEs).  The  open  circles  denote  continuum  FE  nodes  that  have  prescribed  degrees 
of  freedom  (dofs)  D  based  on  the  underlying  grain-scale  response,  while  the  solid  circles 
denote  continuum  FE  nodes  that  have  free  dofs  D  governed  by  the  micromorphic  continuum 
model.  We  intentionally  leave  an  ‘open  window’  (i.e. ,  Direct  Numerical  Simulation  (DNS) 
region)  on  the  particulate  micro-structural  mesh  in  order  to  model  dynamic  failure.  If  the 
continuum  mesh  overlays  the  whole  particulate  micro-structural  region,  as  in  Klein  and  Zim¬ 
merman  [2006]  for  atomistic-continuum  coupling,  then  the  continuum  FEs  would  eventually 
become  too  deformed  by  following  the  micro-structural  motion  during  fragmentation.  The 
blue-dashed  box  at  the  bottom-center  of  the  illustration  is  a  micromorphic  continuum  FE 
region  that  can  be  converted  to  a  DNS  region  for  adaptive  high-fidelity  material  modeling 
as  the  projectile  penetrates  the  target. 


r~\ 


deformable  solid  body  (projectile) 
continuum  FE  mesh 


V  bound  particulate  material  (target) 
multi-scale  computational  model 


-1 

— 

' 

l | 

|  /  f 

particulate  micro-structural  DNS  region  - 

r-'-'7~M==C  /  (pE-and^r'FE/CSEfwY  \  V 

s _ _ /  -tv4-,  / _  _ 1  /M _ A _ --4 — If 

±2. i 

H 

1 1 

1  r 

,  L 

<  r 

to 

coup 
(micromorp 
^articulate  n 

ing  region 
lie  continuum 
iicro-structur< 

FE 

tl  DNS) 

microm 

Drphic  cor 

itinuum  FE 

:  region 

Figure  1.  2D  illustration  of  concurrent  computational  multi-scale  modeling  approach  in  the  contact 
interface  region  between  a  bound  particulate  material  (e.g.,  ceramic  target)  and  deformable  solid 
body  (e.g.,  refractory  metal  projectile). 
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In  order  to  simplify  the  problem,  we  consider  the  compression  and  shear  of  the  interface 
of  DNS  region  with  the  micromorphic  continuum  FE  region  as  a  ID  overlap  problem  between 
glued,  nonlinear  elastic  DEs,  and  micropolar  elastic  continuum  FE  [Regueiro  and  Yan,  2013]. 
Details  may  be  found  in  Regueiro  and  Yan  [2013].  In  summary,  we  learned  that  the  coupling 
methodology  generally  works  for  quasi-statics  between  DEs  and  FEs,  but  the  ID  problem 
is  too  “discrete”  to  yield  the  results  we  expect,  i.e. ,  no  influence  of  overlap  coupling  region. 
This  is  because  shear  and  rotational  degrees  of  freedom  are  coupled  in  the  DE  and  FE 
formulations,  and  thus  upon  shearing  the  end  of  the  DE  string  of  particles,  a  ratcheting 
rotational  mode  is  propagated  through  the  DE  string  and  into  the  FE  continuum  domain 
where  it  is  smoothed  out.  We  are  working  on  extending  this  coupling  method  to  include 
inertia  terms  (dynamics),  and  also  to  couple  grain-FE  to  micromorphic  FE  for  2D  and  3D 
problems. 
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3  Formulation  of  overlap  coupling  for  grain-FE  scale 
and  micromorphic  continuum  FE  scale 

Two  coupling  approaches  in  Sections  3.1  and  3.2  were  considered  in  this  research.  Cou¬ 
pling  Approach  1  is  implemented  with  numerical  results  discussed  in  Section  2.  Coupling 
Approach  1  is  based  on  the  “bridging-scale  decomposition”  in  Wagner  and  Liu  [2003],  Kad- 
owaki  and  Liu  [2004],  Klein  and  Zimmerman  [2006],  and  Coupling  Approach  2  is  based  on 
the  “bridging  domain  method”  in  Xiao  and  Bclytschko  [2004] .  Coupling  Approach  1  couples 
the  DNS  and  micromorphic  FE  dofs  directly  through  interpolation  functions  in  the  overlap 
region,  and  partitions  the  energy  through  scaling  factors  within  each  term  of  the  energy 
(kinetic  and  potential  energies).  Coupling  Approach  2  partitions  the  energies  through  a 
scalar  distance  function  ou(X)  (describing  the  width  of  the  overlapping  domain),  and  uses 
Lagrange  multipliers  to  constrain  the  DNS  displacements  to  the  micromorphic  continuum 
FE  displacements  in  the  overlapping  domain. 

3.1  Coupling  Approach  1: 

It  is  assumed  that  the  micromorphic  continuum-FE  mesh  covers  the  domain  of  the  problem 
in  which  the  bound  particulate  mechanics  is  not  dominant,  whereas  in  regions  of  significant 
grain-matrix  debonding  or  intra-granular  cracking  leading  to  a  macro-crack,  a  grain-scale 
mechanics  representation  is  used  (e.g.,  grain-FE,  grain-DE-FE,  grain  peridynamics  [Silling, 
2000])  in  Fig. 2. 

3.1.1  3D  Kinematics 

Here,  a  summary  of  the  kinematics  of  the  coupled  regions  is  given  for  general  3-D  kinematics, 
following  the  illustration  shown  in  Fig. 2.  It  is  assumed  that  the  micromorphic  continuum 
finite  element  (FE)  mesh  covers  the  domain  of  the  problem  in  which  the  material  is  not 
fragmenting,  whereas  in  regions  of  fracture  and  fragmentation,  a  DNS  (e.g.,  DE/FE/CSE, 
or  peridynamics  [Silling,  2000])  representation  can  be  used.  In  Fig.2,  discrete  domains  are 
defined,  where  the  yellow  background  denotes  the  FE  overlap  region  Bh  with  underlying 
ghost  DNS  nodes,  brown  the  micromorphic  FE  continuum  region  Bh  with  no  underlying 
DNS,  and  white  background  the  free  DNS  region  Bh  UBDNS .  In  summary,  the  finite  element 
domain  Bh  is  the  union  of  pure  micromorphic  continuum  FE  domain  Bh,  overlapping  FE 
domain  with  underlying  ghost  DNS  nodes  Bh,  and  overlapping  FE  domain  with  underlying 
free  DNS  nodes  Bh,  such  that  Bh  =  BhUBhL)Bh .  The  pure  DNS  domain  with  no  overlapping 
FE  domain  (i.e. ,  the  ‘open-window’)  is  indicated  by  BDNS .  The  goal  is  to  have  the  overlap 
region  Bh  U  Bh  as  close  to  the  region  of  interest  (e.g.,  penetrator  skin)  as  to  minimize  the 
number  of  DNS  nodes,  and,  thus,  computational  effort.  Following  some  of  the  same  notation 
presented  in  Klein  and  Zimmerman  [2006],  we  define  a  generalized  dof  vector  Q  for  DNS 
displacements  in  the  system  as 

Q  =  \_Qai  ■  ■  ■  i  <?7]r>  a,P,...,'yeA  (1) 

where  qa  is  the  displacement  vector  of  DNS  node  a,  and  A  is  the  set  of  all  DNS  nodes. 
Likewise,  the  micromorphic  finite  element  nodal  displacements  and  micro-displacements  are 
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Figure  2.  Two-dimensional  illustration  of  the  Coupling  Approach  1  between  grain-FE  DNS  and 
micromorphic  FE  continuum  regions. 


written  as 

D  =  [da,db,  ...,dc,0d,0e,...,0/]T  (2) 

a,b, . . .  ,c  E  J\  ,  d,  e, . . . ,  f  e  M. 

where  da  is  the  displacement  vector  of  node  a ,  <pd  is  the  micro-displacement  tensor  in  vector 
form  of  node  d,  N  is  the  set  of  all  nodes,  and  A4  is  the  set  of  micromorphic  finite  element 
nodes  with  micro-displacement  tensor  degrees  of  freedom,  where  A4  C  JV.  In  order  to  satisfy 
the  boundary  conditions  for  both  regions,  the  motion  of  the  DNS  nodes  in  the  overlap  region 
(referred  to  as  “ghost  DNS  nodes,”  cf.  Fig. 2)  is  prescribed  by  the  micromorphic  continuum 
displacement  and  micro-displacement  fields,  and  written  as 

Q  =  iQa,  Qp,  ■  ■  ■ ,  a,/3,...,7  e  A,  A.  G  Bh  (3) 

while  the  unprescribed  (or  free)  DNS  nodal  displacements  are 

Q=\qs,q„----%}T  ,  A£BhUBDNS  (4) 
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where  A  U  A  =  A  and  A  D  A  =  0.  Likewise,  the  displacements  and  micro-displacement 
tensor  components  of  nodes  overlaying  the  DNS  region  are  prescribed  by  the  DNS  motion 
and  written  as 


D  =  [da,db,...,dc,(f)d,(t>e,...,(j)f]T  (5) 

a,b,  ■  ■  ■  ,c  G  Af  ,  d,e, . . . ,  /  e  M 
Bh  U  Bh 

while  the  nnprescribed  (or  free)  nodal  displacements  and  micro-displacements  are 

D  =  \dm,  dn , .  •  • ,  ds,  <fit,  (j)u,  ■  ■  ■ ,  (6) 

m,n, . . . ,  s  e  Af  ,  t,  u, . . . ,  v  e  A4 

A f,Me  Bh  u  Bh 

where  AfuAf  =  Af,AfnAf  =  0,  A4UA4  =  M,  and  M  <1  M  =  0.  Referring  to  Fig. 2, 
the  prescribed  DNS  nodal  dofs  Q  can  be  viewed  as  boundary  constraints  on  the  free  DNS 
region,  and  likewise  the  prescribed  micromorphic  finite  element  nodal  displacements  and 
micro-displacements  D  can  be  viewed  as  boundary  constraints  on  the  finite  element  mesh 
in  the  overlap  region. 

In  general,  the  displacement  vector  of  a  DNS  node  a  can  be  represented  by  the  fi¬ 
nite  element  interpolation  of  the  micromorphic  continuum  macro-displacement  field  uh  and 
micro- displacement  tensor  field  <&h  evaluated  at  the  DNS  node  in  the  reference  configuration 
X such  that 

*D)  =  £ KPdKW  ,  =  E NPCMW  aeA  (7) 

aGJ\f  b£M 

where  N/  are  the  shape  functions  associated  with  the  micromorphic  continuum  displacement 
held  uh,  Nb  are  the  shape  functions  associated  with  the  micromorphic  continuum  micro¬ 
displacement  held  and  h  indicates  the  characteristic  length  of  the  micromorphic  FE 
mesh.  Recall  that  N/  and  N/  have  compact  support  and,  thus,  are  only  evaluated  for  DNS 
nodes  that  lie  within  a  micromorphic  FE  containing  nodes  a  and  b  in  its  domain. 

We  now  derive  the  micro-continuum  displacement  (u')h  of  a  micro-element  (Fig.3)  using 
the  interpolated  macro-displacement  uh  and  micro-displacement  tensor  <l?h  in  Eq.(7).  The 
position  of  DNS  node  a  in  the  current  configuration  is 

(x')h(Xha,  Eha,  t )  =  xh(XhQ,  t )  +  t(Xhal  E*,  t)  (8) 

where  is  the  reference  relative  position  vector  of  DNS  node  a  with  respect  to  reference 
continuum  position  Xha  (e.g.,  a  Gauss  point)  in  £>q,  and  $,h  is  the  current  relative  position 
vector  of  DNS  node  a  with  respect  to  the  current  continuum  position  xh  in  Bh .  Introducing 
the  relations, 


*  {X a,  *) 

(Xha,Bha,t) 


Xha  +  uh(Xha,t) 

Xh(Xha,t)-Zha  =  (l  +  *h(Xha,t))-Zha 


and  substituting  into  Eq.(8),  we  have  the  current  position  of  DNS  node  a  as 


(9) 

(10) 


Final  Report  for  W911NF-09-1-0111:  period  05/01/2009  to  09/30/2012 


7 


Figure  3.  Multiplicative  decomposition  of  deformation  gradient  F  and  micro-deformation  tensor 
X •  Geometrical  points  (“macro-elements”)  with  centroids  C,  C ,  and  c  live  in  their  respective 
configurations:  P  £  £>o,  P  £  B,  and  p  £  B.  Material  points  (“micro-elements”)  with  centroids 
C',  C',  and  d .  Differential  vectors  and  deformable  directors  are  mapped  accordingly:  dx  =  FdX, 
dx  =  FedX ,  dX  =  FpdX ,  £  =  *3,  £  =  xeS,  and  5  =  XPS. 

(x')h(K,  SS,  *)  =  X«  +  t)  +  (1  +  t))  ■  Sha  (11) 

If  we  rewrite  (x')h  in  terms  of  its  reference  position  (. X')h ,  and  we  introduce  a  micro¬ 
continuum  displacement  vector  ( u')h ,  we  have 

(*T t)  =  (X')\Xha,  Sq)  +  (u')h  (X^,  Sha,  t)  (12) 

where  (X')h(^,  Eha)  =  Xha  +  Eha,  then  we  have 

{u')h{Xhal  Ej,  t)  =  (x')h(Xha,  S",  t)-Xha-  Eha  (13) 

Finally,  we  have  the  micro-continuum  displacement  vector  ( u')h  associated  with  DNS  node 
a,  such  that 

qa(t)  =  e u')h(Xha,~ha,t )  =  uh(Xha,t )  +  4>h(-X*t)  •  (14) 

Then,  for  all  ghost  DNS  nodes  (cf.  Fig. 2),  the  interpolations  can  be  written  as 

Q  =  XqD  ■  D  +  ■  D  (15) 

where  N qD  and  AT" gg  are  shape  function  matrices  containing  individual  nodal  shape  func¬ 
tions  N/  and  IV*,  but  for  now  these  matrices  will  be  kept  general  to  increase  our  flexibility 
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in  choosing  interpolation/projection  functions  (such  as  those  used  in  meshfree  methods). 
Overall,  the  DNS  nodal  displacements  may  be  written  as 

Q 

Q 

where  Q'  is  introduced  [Klein  and  Zimmerman,  2006]  as  the  error  (or  “fine-scale”  [Wagner 
and  Liu,  2003])  in  the  interpolation  of  the  free  DNS  nodal  displacements  Q,  whose  function 
space  is  not  rich  enough  to  represent  the  true  free  DNS  motion.  The  shape  function  matrices 
N  are,  in  general,  not  square  because  the  number  of  free  DNS  dofs  are  not  the  same  as  free 
and  prescribed  micromorphic  FE  dofs,  and  the  number  of  ghost  DNS  dofs  is  not  the  same 
as  prescribed  and  free  micromorphic  FE  dofs.  A  scalar  measure  of  error  e  in  DNS  nodal 
displacements  is  defined  as  [Klein  and  Zimmerman,  2006] 


N 


QD 


N 


QD 


n6d  n. 


D 

+ 

Q' 

D 

0 

(16) 


e  =  Q'  ■  Q’ 


(17) 


which  may  be  minimized  with  respect  to  prescribed  micromorphic  continuum  FE  nodal  dofs 
D  to  solve  for  D  in  terms  of  free  DNS  dofs  and  micromorphic  continuum  FE  nodal  dofs  as 

5  =  MBBNId(Q  -  NQ°D)  .  MdB  =  NtqBNqB  (18) 

This  is  known  as  the  “discretized  L2  projection”  [Klein  and  Zimmerman,  2006]  of  the  free 
DNS  nodal  dofs  Q  and  free  micromorphic  FE  nodal  dofs  D  onto  the  prescribed  micromorphic 
FE  nodal  dofs  D.  Upon  substituting  Eq.(18)  into  Eq.(15),  we  may  write  the  prescribed 
DNS  dofs  Q  in  terms  of  free  DNS  dofs  Q  and  micromorphic  continuum  FE  nodal  dofs  D. 
In  summary,  these  relations  are  written  as 


where 


Q 

—  bqqQ  +  BQDB 

(19) 

D 

=  BdqQ  A  Bq dD 

(20) 

Bqq 

=  NqqBqq 

(21) 

bqd 

—  Nqd  +  NqqBqd 

(22) 

Bdq 

=  MAaA 

DD  QD 

(23) 

Bdd 

=  -M^N^NqD 

DD  QD  ^ 

(24) 

As  shown  in  Fig. 2,  for  a  finite  element  implementation  of  this  dof  coupling,  we  expect  that 
free  DNS  dofs  Q  will  not  fall  within  the  support  of  free  micromorphic  continuum  FE  nodal 


dofs  D ,  such  that  it 

can  be  assumed  that  Nqd  =  0  and 

Q  =  B^qQ  +  B^dD  ,  D  =  BqqQ 

(25) 

where 

BqQ  =  NqQBQq  ,  Bqd  =  iVgD 

(26) 

Bdq  =  Md1dNQD  ’  Bdd  =  0 

(27) 

The  assumption  Nqd  /  0  would  be  valid  for  a  meshfree  projection  of  the  DNS  nodal  dofs 
to  the  micromorphic  FE  nodal  dofs,  as  in  Klein  and  Zimmerman  [2006],  where  we  could 
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imagine  that  the  domain  of  influence  of  the  meshfree  projection  could  encompass  a  free 
DNS  node;  the  degree  of  encompassment  would  be  controlled  by  the  chosen  support  size  of 
the  meshfree  kernel  function.  The  choice  of  meshfree  projection  in  Klein  and  Zimmerman 
[2006]  was  not  necessarily  to  allow  Q  be  projected  to  D  (and  vice  versa),  but  to  remove  the 
computationally  costly  calculation  of  the  inverse  in  Eqs.  (19)  and  (20). 

3.1.2  3D  Kinetic  and  Potential  Energy  Partitioning  and  Coupling 

For  the  DNS  equations,  the  kinetic  energy  is  T®,  dissipation  function  F® ,  and  potential 
energy  U®,  such  that 


Tq  =  i QMqQ 

FQ  =  aQTQ  (28) 

rQ 

Uq(Q)  =  /  F1NT,Q(S)dS 

Jo 

where  is  the  mass  matrix,  and  FINT,<^  is  an  internal  force  vector  for  the  DNS  region. 
The  dissipation  function  F ®  is  written  as  a  linear  function  of  the  kinetic  energy  T ®  with 
proportionality  coefficient  czq,  which  falls  within  the  class  of  damping  called  “Rayleigh  damp¬ 
ing”  (pg.130  of  Rayleigh  [1945]).  For  the  micromorphic  continuum  FE  equations,  TD  is  the 
kinetic  energy,  FD  the  dissipation  function,  and  UD  the  potential  energy,  such  that 


rpD  _ 


Fd  = 


-DMdD 

2 

aDTD 

,D 


Ud(D) 


?INT,D 


1 0 


(S)dS 


(29) 


where  M°  is  the  mass  and  micro-inertia  matrix,  and  F 1 NT'  °  is  the  internal  force  vec¬ 
tor.  We  assume  the  total  kinetic  and  potential  energy  and  dissipation  of  the  coupled  DNS- 
micromorphic-continuum  system  may  be  written  as  the  sum  of  the  energies 


T(Q,D) 

=  T«(Q,Q(Q,D))  +  Td(D,D(Q)) 

(30) 

U(Q.D ) 

=  UQ(Q,Q(QjD))  +  Ud(D.D(Q)) 

(31) 

F(Q.D) 

=  Fq(Q.^(Q,D))  +  Fd(D,D(Q)) 

(32) 

where  we  have  indicated  the  functional  dependence  of  the  prescribed  DNS  dofs  and  mi¬ 
cromorphic  FE  dofs  solely  upon  the  free  DNS  dofs  and  micromorphic  FE  dofs  Q  and  D. 
respectively.  Lagrange’s  equations  may  then  be  stated  as 


d  fdT\  d T  OF  dU 
dt  \dQj  ~  dQ  +  dQ  +  dQ 
d  ( dT\  dT  dF  dU 
dt\dDj  ao  +  ao  +  ao 


piEXT.Q 

p*EXT,D 


(33) 
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which  lead  to  a  coupled  system  of  governing  equations  (linear  and  first  moment  of  momen¬ 
tum)  for  the  coupled  DNS-micromorphic-continuum  mechanics.  The  derivatives  are 


dT 

dQ 

dF 

OQ 

dU 

dQ 

dT 

~dD 

dF 

dt> 

dU 

dD 


dTQ  dTQ 


dTD 


„  dT  n 
’  dQ~° 


an 


=  0 


std 

dD 


B 


DQ 


fSTD\ 

\~8d) 


(34) 

(35) 

(36) 

(37) 

(38) 

(39) 


If  the  potential  energy  U  is  nonlinear  with  regard  to  DNS  constitutive  models  and  rnicro- 
morphic  elasto-plasticity,  then  Eq.(33)  may  be  integrated  in  time  and  linearized  for  solution 
by  the  Newton- Raphson  method.  The  benefit  of  this  multiscale  method,  as  pointed  out  by 
Wagner  and  Liu  [2003],  is  that  time  steps  may  be  different  for  the  DNS  and  micromorphic 
FE  solutions.  A  multiscale  time  stepping  scheme  will  follow  an  approach  similar  to  Wagner 
and  Liu  [2003].  To  complete  Eq.(33)  and  identify  an  approach  to  energy  partitioning,  the 
individual  derivatives  may  be  written  as 


dTQ 

~dQ 

dTD 

dD 

dF 

dQ 

dF 

~&D 

dUQ 

~dQ 

dUD 

dD 


ary  -  a 

=  MqQ  ,  — v-  =  MqQ 
dQ 

dT°  ~  a 

=  m'  d  ,  — =  MdD 

dD 

aQ  (M«Q  +  +  aD  (. BTBQMBb ) 


aQ(BTdnMQQ)+aD  (Md D 
dUQ 


QD 

Fint’q(Q )  , 


dQ 

dUD 


_  jpINT,Q 


(Q) 


MINT’D (£))  ;  =  FINT’°(D) 

db 


(40) 

(41) 

(42) 

(43) 

(44) 

(45) 


where  superscript  Q  denotes  free  DNS  dofs  and  Q  ghost  DNS  dofs,  whereas  superscript  D 
denotes  free  micromorphic  FE  dofs  and  D  prescribed  micromorphic  FE  dofs.  The  energy 
partitioning  will  be  introduced  through  the  definition  of  these  terms  below.  First,  substitute 
Eqs.  (40)-(45)  into  Eq.(33)  to  arrive  at  the  coupled  nonlinear  equations  in  terms  of  Q  and 
D  as 
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(m«  +  BtQqM$B3q  +  BlQMBBBQ )  Q 
+  {bTqQmS>Bqd  +  BtBqMBBBd)  D 

+  (®qMq  +  aQB3QMSB3q  +  apB  BqMB  B/  Q  +  uqB qQMB B 3dD 
+F,nt*(Q)  +  BtSqF'nt-B{B3qQ  +  B3dD] 

+Bt3qFint’3[B3qQ  +  BBdD }  =  Fext’q  +  B/fQFEXT&  (46) 

(■ B8dmSB3d  +  BtBdMBBBq)  Q 
+  (md  +  bt3dmbb3d  +  btBdmbbBd)  D 

+aQB 3dMB B 3qQ  +  (aQB3DMBB3D  +  aBM  )  D 

+bt3df,nt,bib3qq  +  BScD]  +  BtBdF,nt'd[BB3Q  +  BBdD] 

+F,nt’d(D)  =  Fext-d  +  btBdFextB  (47) 

where  an  expression  in  brackets  [•]  denotes  the  functional  dependence  of  the  nonlinear  inter¬ 
nal  force  vector.  Note  the  projections  through  the  B  matrices  of  the  corresponding  mass  and 
damping  matrices,  and  forcing  vectors.  First,  starting  with  the  mass  matrices  for  the  DNS 
nodes,  we  partition  the  kinetic  energy  (and,  in  turn,  the  dissipation  functions)  as  follows: 

M3  =  (1  -  f )  m®  ,  aei,  xa  G  Bh  (48) 

OL 

Mq  =  Mdns’q  +  MQ  (49) 

mdns,q  =  fa  m/j  ,  p  e  A  ,  Xp  e  BDNS 

MQ  =  (1  —  f)  ,  P  e  A  ,  xp  e  Bh 

p 

where  is  the  mass  matrix  of  ghost  DNS  nodes  in  Bh,  MDNS the  mass  matrix  of  free 

DNS  nodes  in  BDNS ,  M  the  mass  matrix  of  free  DNS  nodes  in  Bh,  and  f  is  a  weighting 
factor  for  the  kinetic  energy  in  the  overlap  region  Bh  U  Bh .  For  no  homogenized  micromor- 
phic  continuum  contribution  to  the  kinetic  energy  in  the  overlap  region,  f  =  0,  and  for  full 
micromorphic  continuum  homogenization  of  the  underlying  DNS  kinetic  energy,  f  —  1.  In 
our  case,  we  will  consider  the  range  0  <  f  <  1.  Given  that  the  proposed  multiscale  modeling 
framework  is  to  be  used  in  an  adaptive  fashion  in  the  future,  having  an  overlaying  micromor¬ 
phic  continuum  homogenization  of  the  DNS  response  is  attractive  when  DNS  is  converted  to 
micromorphic  continuum  representation,  and  vice  versa  (in  a  statistical  manner*).  For  the 

*  statistical,  in  the  sense  of  generating  a  DNS  representation  from  a  micromorphic  continuum  one,  where 
the  underlying  DNS  region  does  not  exist;  converting  from  DNS  to  micromorphic  continuum  representation 
is  straightforward  given  the  built-in  homogenization  that  the  micromorphic  continuum  possesses 
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mass  and  micro-inertia  matrices  associated  with  the  micromorphic  continuum,  we  partition 
the  kinetic  energy  as  follows: 


Md  = 


—  D  - — D 

M  +  M 


(50) 


—  D 

M  = 


M  0 


0  M 

M"  =  J\(f  (mu,e)  +  femu’e) 


e&Bh 


M  —  (f  +  rem$,e) 


eg  Bh 


- — D 

M  =  f 


M  0 


0  M 


Md  = 


' — '  D  —  n 

M  +M 


MD  = 


M  0 
0  M* 


,  MU  =  ^  (mu,e)  ,  M  =  A  {m*'e) 

eg  Bh  eg  Bh 


,  M"=^  mue  ,  =  A  m$,e 

e£@h  e€@h 


(51) 


where  M  is  the  micromorphic  continuum  mass  and  micro-inertia  matrix  associated  with 

~  —  D 

prescribed  micromorphic  FE  dofs  in  B'.  M  the  micromorphic  continuum  mass  and  rnicro- 

—  5  - — D 

inertia  matrix  associated  with  free  micromorphic  FE  dofs  in  £>  ,  where  M  and  M  are 

extracted  from  the  total  mass  and  micro- inertia  matrix  M  in  Bh ,  with  superscript  (•)^> 
denoting  the  full  mass  and  micro-inertia  matrix  associated  with  elements  in  Bh,  (•)  is  a 
homogenization  operator,  fe  is  the  partitioning  coefficient  of  micromorphic  continuum  kinetic 
energy  associated  with  element  Be  C  Bh.  A  simple  choice  is  a  volume  fraction  fe  =  Be,D /Be, 
where  Be,D  =  Be  —  Be’®]  B/,D  is  the  non-overlapping  micromorphic  continuum  part  of  element 
volume  Be  C  Bh,  and  Be’®  is  the  overlapped  ghost  DNS  volume  in  the  element  (cf.  Fig. 2).  For 
kinetic  energy  partitioning,  a  volume  fraction  that  directly  relates  to  mass  and  micro- inertia 

. —  D 

partitioning  seems  an  appropriate  choice.  M  is  the  homogenized  micromorphic  continuum 
mass  and  micro-inertia  matrix  associated  with  prescribed  micromorphic  FE  dofs  in  Bh\  where 
if  f  =  0,  there  is  no  micromorphic  continuum  homogenization  in  Bh  (i.e. ,  all  kinetic  energy  is 
due  to  underlying  DNS).  MD  is  the  micromorphic  continuum  mass  and  micro-inertia  matrix 
associated  with  free  micromorphic  FE  nodal  dofs  in  the  pure  micromorphic  continuum  FE 
domain  Bh. 

For  the  potential  energy  (internal  force)  and  external  force  partitioning  in  the  DNS 
system,  we  write 
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jpINT,Q  _ 

(1  -  q)  A  flNT,Q  ,xeeBh 

(52) 

p!NT,Q  _ 

jpINT,DNS,Q  _|_  fINT’Q 

(53) 

FINT,DNS,Q  =  A  fINT,Q  xs  e  BDNS 
<5 


FINT,Q  =  (1  -  q)  A  /f  T’Q  ,  G  B 


<5 

fEXT,Q 

=  (l~q)  /\f?XT’Q  ,  XetBh 

(54) 

fEXT,Q 

6 

j?EXT,DNS,Q  fEXT’Q 

(55) 

FEXT,DNS,Q  =  ^  fEXT,Q  ^  £  jgDNS 
£ 

FEXT'Q  =  (1  -  q)  A  f!XT’Q  ,  Xs  e  Bh 

<5 


where  F1NT’Q  is  the  internal  force  vector  associated  with  ghost  DNS  nodes  in  Bh  ,  j?INT-DNS& 
is  the  internal  force  vector  associated  with  free  DNS  nodes  in  BDNS ,  F  is  the  internal 
force  vector  associated  with  free  DNS  nodes  in  Bn,  F  ,w  is  the  external  force  vector  asso- 

~  j  77 1  ~\/'rn  n  n  /o 

ciated  with  ghost  DNS  nodes  in  Bh  ,  F  ’  is  the  external  force  vector  associated  with 
free  DNS  nodes  in  BUNb ,  F  is  the  external  force  vector  associated  with  free  DNS  nodes 
in  Bh,  and  q  is  a  weighting  factor  for  the  potential  energy  in  the  overlap  region  Bh  U  Bh .  For 
no  homogenized  micromorphic  continuum  contribution  to  the  potential  energy  in  the  overlap 
region,  q  —  0,  and  for  full  micromorphic  continuum  homogenization  of  the  underlying  DNS 
potential  energy,  q  —  1.  In  our  case,  we  will  consider  the  range  0  <  q  <  1.  Note  that  in  Klein 
and  Zimmerman  [2006],  they  chose  q  =  0.  Their  Cauchy-Born  elastic  constitutive  model  acts 
like  a  homogenization  operator  on  the  underlying  atomistic  response,  but  instead  of  keeping 
an  overlain  Cauchy-Born  representation,  the  potential  energy  is  completely  represented  by 
the  underlying  atomistic  response,  except  in  the  overlap  region  Bh  where  partitioning  occurs. 
Note  that  xe  and  x$  denote  positions  of  DNS  nodes  for  calculating  internal  force  vectors  in 
Eqs.  (52)-(56),  whereas  xa  and  xg  in  Eqs.  (48)  and  (49)  denote  DNS  nodes  for  calculating 
DNS  mass  matrices. 

For  the  potential  energy  (internal  force)  partitioning  in  the  micromorphic  continuum,  we 
write 

(56) 


j;INT,u,e  pINT,<S>  _  j, INT,<S>,e 

e£Bh  eG@h 


jpINT,D  _  fINT’D  pINT,D 


INT,D 


rINT,u 


and 
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jpINT,D 


~  INT,D  —  INT,D 

F  +F 


(57) 


- 1  NT,D 


-  INT,u 


•  lNT,u 


-  I  NT,® 


A  0(f 

e&Bh 


INT,u,e 


>  +  ?fINT™) 


INT,u,e\  _ 


f 


■  I  NT,® 


)=  I  (B 


u,e\T 


A  (Hf 

e£Bh 

INT,u 


<T )  dv 

/jVT,4>,e^  _|_  ~e 


^  I  NT, D 

F  =  q 


- 1  NT,u 


A 

e&Bh 


—  INT,® 

F 

^INT,u,e\^ 


INT,$ 


A  </ 

e&Bh 


INT,®,e\ 


~INT,D 

where  is  the  internal  force  vector  associated  with  free  micro morphic  FE  dofs  in  Btl , 

~  INT,D 

E1  the  internal  force  vector  associated  with  prescribed  micromorphic  FE  dofs  in 

~  I  NT, D  ~INT,D  ~  INT,D 

where  -F  and  F  are  extracted  from  the  full  internal  force  vector  F  ,  with 
superscript  (•)-D  denoting  the  full  internal  force  vector  associated  with  elements  in  Bh,  qe 
is  the  partitioning  coefficient  of  micromorphic  continuum  potential  energy  associated  with 
element  Be  C  Bh,  and  (•)  is  a  homogenization  operator  (to  be  defined  later).  A  simple  choice 
is  a  volume  fraction  qe  =  fe.  Klein  and  Zimmerman  [2006]  considered  a  more  sophisticated 
approach  using  an  atomic  bond  density  function  solved  to  reproduce  a  minimum  potential 
energy  state  for  homogeneous  deformation.  The  analogy  here  for  DNS  nodes  would  be  a  DNS 
potential  energy  term  (internal  force  vectors).  This  will  be  considered  further  in  future  work. 
For  now,  we  consider  a  volume  fraction  partitioning  through  qe ,  and  a  simple  scaling  through 

—  INT,D 

coefficients  q.  F  is  the  homogenized  internal  force  vector  associated  with  prescribed 
micromorphic  FE  dofs  in  Bh ,  which  has  no  contribution  if  q  —  0,  i.e.,  underlying  DNS  internal 
forces  provide  full  contribution  in  Bh.  F  '  is  the  internal  force  vector  associated  with 
free  micromorphic  FE  dofs  in  the  pure  micromorphic  continuum  domain  Bh .  The  external 
force  vector  is  written  as 


n  EXT,D 


~EXT,D  pEXT,D 


-,EXT,D 


EXT,u 


-.EXT,® 


pEXT,u  =  Ft  +  Fu  +  ^  f 

e£@h 


EXT,u,e 

b 


fext,®  = 


A  f 

e£Bh 


EXT,®,e 

e 


(58) 
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and 


jpEXT.D 


~EXT,D  — EXT,D 

F  +  F 


-  EXT,D 


-  EXT,u 


-  EXT,® 


■  EXT,u 


-  EXT,® 


A  [4{fb 

e£Bh 


EXT,u,e 


i  ~e  nEXT.u.e 

+ <i  /., 


A 

e£Bh 


q(f 


EXT,®,e 


+  qef 


e  _pEXT,&,e 


^ EXT,D 

F  —  q 


■  EXT,u 


■  EXT,® 


■  EXT,u 


A  </ 

e£Bh 


EXT,u,e 

b 


■EXT,® 


A 

e£Bh 


(59) 


/ 


EXT,®,e 

l 


~  EXT,D 

where  -F  is  the  external  body  force  and  couple  vector  associated  with  free  micromorphic 

~EXT,D 

FE  nodal  dofs  in  F  the  external  body  force  and  couple  vector  associated  with 

~  EX T,D  ~  EXT,D 

prescribed  micromorphic  FE  dofs  in  B  ,  where  F  and  F  are  extracted  from 

~  EXT,D  ^  EXT,D 

F  ,  the  total  external  body  force  and  couple  vector  calculated  in  Bh.  F  is  the 

homogenized  external  body  force  and  couple  vector  associated  with  prescribed  micromorphic 
FE  dofs  in  Bh,  which  has  no  contribution  if  q  —  0,  i.e.,  underlying  DNS  body  forces  provide 
full  contribution  in  Bh .  F  is  the  external  force  and  couple  vector  associated  with  free 
micromorphic  FE  dofs  in  the  pure  micromorphic  continuum  FE  domain  Bh . 


3.2  Coupling  approach  2: 

In  this  coupling  approach,  it  is  also  assumed  that  the  micromorphic  continuum-FE  mesh 
covers  the  domain  of  the  problem  in  which  the  bound  particulate  mechanics  is  not  dominant, 
whereas  in  regions  of  significant  grain-matrix  debonding  or  intra-granular  cracking  leading 
to  a  macro-crack,  a  grain-scale  mechanics  representation  is  used  in  Fig.4. 

In  Fig.4,  discrete  domains  are  defined,  where  the  yellow  background  denotes  the  DNS- 
FE  overlap  region  ffDNS-FE ^  with  DNS  region  denoted  by  BDNS  terminating  at  dBDNS ,  and 
micromorphic  continuum  FE  region  denoted  by  BhE  terminating  at  dBh  E . 

3.2.1  3D  Kinematics 

The  notation  here  for  Coupling  Approach  2  is  similar  to  that  used  for  Coupling  Approach 
1.  We  define  a  generalized  dof  vector  Q  for  all  DNS  nodal  displacements  in  the  system  as 

Q  =  \q.oti  Qf 3,  •  ■  • ,  <?7]T,  a,P,...,yeA  (60) 

where  qa  is  the  displacement  vector  of  DNS  node  a,  and  A  is  the  set  of  all  DNS  nodes. 
Likewise,  the  micromorphic  finite  element  nodal  displacements  and  micro-displacements  are 
written  as 
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bound 

particulate 

material 

target 


DNS  (DE/FE/CSE) 


coupling 
between  DNS 
and 

micromorphic 
continuum  FE 


micromorphic 
continuum  FE 


shear/compression  loading 


bdns 

\K-dB™ 

\bdns-fe 

1*.  Jt  (JB/)XS 
Lfe 


Figure  4.  Two-dimensional  illustration  of  Coupling  Approach  2  between  grain-FE  DNS  and 
micromorphic  FE  continuum  regions. 


D  =  [da,db,...,dc,<i>d,<t)e,...,<t>f]T  (61) 

a,  b, . . . ,  c  E  Af  ,  d,e, . . . ,  /  6  M 

where  da  is  the  displacement  vector  of  node  a,  <pd  is  the  micro-displacement  tensor  in  vector 
form  of  node  d,  A f  is  the  set  of  all  nodes,  and  A4  is  the  set  of  micromorphic  finite  element 
nodes  with  micro-displacement  tensor  degrees  of  freedom,  where  A4  C  Af. 

3.2.2  3D  Kinetic  and  Potential  Energy  Partitioning  and  Coupling 

For  the  DNS  equations,  the  kinetic  energy  is  T®,  dissipation  function  F® ,  and  potential 
energy  UQ ,  such  that 


Tq  =  i QMqQ 

Fq  =  aQTQ  (62) 

rQ 

Uq(Q)  =  /  F1NT,Q(S)dS 

Jo 

where  is  the  mass  matrix,  and  FINT’®  is  an  internal  force  vector  for  the  DNS  region. 
For  the  micromorphic  continuum  FE  equations,  TD  is  the  kinetic  energy,  FD  the  dissipation 
function,  and  UD  the  potential  energy,  such  that 
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Td  =  -DMdD 

2 

F°  =  cldTd  (63) 

rD 

Ud(D)  =  /  FINT’D(S)dS 

Jo 

where  MD  is  the  mass  and  micro-inertia  matrix,  and  FIN1-D  is  the  internal  force  vec¬ 
tor.  We  assume  the  total  kinetic  and  potential  energy  and  dissipation  of  the  coupled 
DNS-micromorphic-continuum  system  may  be  written  as  the  sum  of  the  energies  using  the 
scalar  distance  function  u>(X)  =  f(X)/£0,  where  f(X)  is  the  orthogonal  distance  of  point 
X  E  Bdns~fe  from  dBFE ,  and  f0  is  the  orthogonal  length  of  fiDNS-FE,  such  that 


T(Q.D) 

=  (1  -lu)Tq(Q)+ujTd(D) 

(64) 

U(Q,D) 

=  (1  -oj)Uq{Q)+ojUd(D) 

(65) 

F(Q,D) 

=  (1  -uj)Fq(Q)+uFd{D) 

(66) 

Furthermore,  the  potential  energy  U  of  the  coupled  system  is  modified  to  include  a  constraint 
using  a  Lagrange  multiplier,  such  that  the  Augmented- Lagrange  multiplier  constrained  po¬ 
tential  energy  U\ ugLag  is  written  as 

LAugLag  =  U  +  A  ■  g  +  -kg  ■  g  (67) 

where  A  is  the  vector  of  Lagrange  multipliers  (forces  acting  on  DNS  nodes  in  the  overlap 
region),  g  is  the  vector  of  constraints  applied  at  each  DNS  node  a  within  the  overlap  region, 
and  k  is  a  penalty  parameter.  The  constraint  at  DNS  node  a  in  the  overlap  region  is  written 
as 

9a  =  ( K)h(t )  -  9a(t)  =  0  (68) 

where  ( u'Q)h{t )  is  the  micro-continuum  displacement  at  node  a  shown  in  Eq.(14).  Lagrange’s 
equations  may  then  be  stated  as 

d_  (dT\  _  dT_  dF/  dT/AugLag  _  _  jpEXT,Q 

dt  \dQ )  dQ  OQ  dQ 

d_  (dT_\  _  &r_  dF/  dUA ugLag  _ 

zpEXT,D 

It[dD)-dD  +  ^D  +  -dD-=F  m 

which  lead  to  a  coupled  system  of  governing  equations  (linear  and  first  moment  of  momen¬ 
tum)  for  the  coupled  DNS-micromorphic-continuum  mechanics.  After  taking  the  derivatives, 
and  assuming  mass  and  micro-inertia  lumped  matrices  (for  explicit  dynamics),  we  arrive  at 
the  coupled  system  of  equations  to  solve  at  DNS  node  a 

A  {(1  -  “(VJ)  +  flNT'Q(Q)\  -K-  kga  =  /fx™}  (70) 

a 

where  m®  is  the  lumped  mass  at  DNS  node  a,  fIaNT,(^  is  the  internal  force  vector  at  DNS 
node  a,  f^XT,(^  is  the  external  force  vector  at  DNS  node  a.  For  micromorphic  FE  node  a, 
we  have 
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A 


m°da 


+  aDm f  da  +  / 


INT,D , 
a 


D) 


+  A  K(x«)  (A«  +  *»«)]  =  f 


EXT,D 
a 


(71) 

where  m1/  is  the  lumped  mass  at  micromorphic  FE  node  a,  fIaNT’D  is  the  internal  force  vector 
at  micromorphic  FE  node  a,  faXT’D  is  the  external  force  vector  at  micromorphic  FE  node 
a.  For  micromorphic  FE  node  d,  we  have 


u 


:xn 


id  0 d  +  aDid  4>d  +  f 


:D  . 


fINT,D, 
d  ' 


D) 


/\[N*(Xa)(\a  +  kga)®Zha\  = 


jEXT,D 


(72) 


where  if/  is  the  lumped  micro-inertia  at  micromorphic  FE  node  d. 

An  explicit  time  integration  procedure  and  solution  method  is  described  in  Xiao  and 
Bclytschko  [2004],  which  can  also  be  generalized  for  implicit  time  integration  and  a  Newton- 
Raphson  nonlinear  solution  procedure. 
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4  Overview  of  adaptivity  approach 

Numerical  adaptivity  is  introduced  to  track,  at  the  grain-scale,  regions  of  dynamic  material 
failure.  An  adaptive  scheme  can  be  considered  in  two  ways  to  calculate  conversion  criteria: 
(a)  using  a  calculated  shear/compaction  strain  energy  density  averaged  over  a  RVE  of  DNS 
micro-structure,  and  calculated  at  a  Gauss  point  in  a  micromorphic  FE;  and  (b)  a  posteriori 
model  error  estimators  (e.g.,  Oden  et  al.  [2006]).  In  terms  of  method  development  for 
conversion,  we  consider  the  cases  illustrated  in  Fig. 5  where  a  continuum  micromorphic  FE 
mesh  region  is  converted  to  a  particulate  micro-structural  DNS  region,  and  vice  versa.  One 
of  the  challenges  with  this  conversion  is  that  we  do  not  know  the  material  micro-structure 
to  make  up  the  DNS  region,  nor  its  material  parameters.  We  can  use  probability  density 
functions  (PDFs)  to  sample  approximate  grain  size,  shape,  percentage  of  matrix  versus  grain 
material,  inter-/intra-granular  constitutive  parameters,  and  grain- matrix  interface  de-bond 
strength,  for  instance.  The  sampling  of  PDFs  will  provide  an  initial  guess  for  the  DNS 
micro-structural  geometry  and  material  parameters.  To  ensure  smooth  conversion  in  terms 
of  average  stress,  strain,  and  energy,  we  can  calibrate  inversely  the  material  parameters  (using 
the  PDF  sampling  as  an  initial  guess,  and  the  DNS  geometry  fixed)  through  a  micromorphic 
homogenization  procedure. 


Figure  5.  Adaptivity:  conversion  from  micromorphic  continuum  FE  to  particulate  DNS  material 
model  representation,  and  vice  versa.  The  micromorphic  continuum  FE  region  could  be  the  blue- 
dashed  box  at  the  bottom-center  of  the  illustration  in  Fig.f  to  be  converted  to  a  particulate  DNS 
region  for  adaptive  high-fidelity  material  modeling  as  the  projectile  penetrates  the  target.  Note  that 
since  we  do  not  know  the  exact  particulate  DNS  configuration  and  material  properties  when  con¬ 
verting  from  continuum  FE  to  DNS  in  a  region  of  potential  fracture  nucleation  and  fragmentation, 
this  process  is  probabilistic,  and  thus  PDFs  will  be  used  for  the  conversion. 

We  ended  up  not  having  time  to  implement  these  conversions  illustrated  in  Fig. 5. 
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5  Future  Research 

Future  research  entails  the  following:  (1)  extend  the  micromorphic  3D  FE  implementation 
to  explicit  dynamics  in  Tahoe;  (2)  implement  the  overlap  procedures  presented  in  Sect. 3; 
and  (3)  test  the  adaptivity  concepts  presented  in  Sect. 4. 

Overall,  the  full  3D  finite  strain  micromorphic  coupled  FE  implementation  and  Drucker- 
Prager  model  formulation  and  implementation  took  longer  than  expected,  but  now  it  is 
nearly  completed  for  the  quasi-static  case.  The  overlap  coupling  in  3D  will  continue  after 
this  project,  to  be  supported  by  another  project  and  to  include  the  multiphase  (solid,  liquid, 
gas)  nature  of  geologic  and  particulate  materials.  The  adaptivity  of  our  multiscale  modeling 
approach  will  be  investigated  as  part  of  future  work,  and  it  is  still  an  important  problem  to 
investigate. 
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